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ABSTRACT 


Transport properties of dense liqnid helinm nnder the conditions of planet’s 
core and cool atmosphere of white dwarfs have been investigated by nsing the 
improved centroid path-integral simnlations combined with density fnnctional 
theory. The self-diffnsion is largely higher and the shear viscosity is notably 
lower predicted with the qnantnm mechanical description of the nnclear motion 
compared with the description by Newton eqnation. The resnlts show that nn¬ 
clear qnantnm effects (NQEs), which depends on the temperatnre and density of 
the matter via the thermal de Broglie wavelength and the ionization of electrons, 
are essential for the transport properties of dense liqnid helinm at certain astro- 
physical conditions. The Stokes-Einstein relation between diffnsion and viscosity 
in strongly conpled regime is also examined to display the inflnences of NQEs. 


Subject headings: dense matter - dense helium - diffusion - viscosity - quantum effects 


1. Introduction 


Helium, the second most abundant element in the universe, is one of the main 
components of giant planets and a large number of recently discovered exoplanets (Santos 
et ah 2005; McMahon et ah 2012). Meanwhile, helium is also expected to be largely present 
in the cool white dwarf atmospheres (Winisdoerffer & Chabrier 2005), which generally 
affect the cooling rate of white dwarfs (WDs) when having exhausted their nuclear fuel. 
Matter in these astrophysical bodies is generally under extreme temperature and pressure 
conditions (Chabrier et ah 2007). For instance, the temperature range of Jupiter, a typical 
solar giant planet, is estimated from 6500 K to 20000 K, while the pressure can reach 70 
Mbar (Fortov 2011). For some exoplanets, it can be up to 19 Gbar (Swift et ah 2012). 
The density range of the outer layer of WDs is a few g/cm^ and the temperature retains 
thousands of Kelvin (Iglesias et al. 2002). The central element for determining the structure 
and evolution of these astrophysical objects is the accurate description of the physical 
properties of matter under the relevant conditions of temperature and pressure (Fortney 
& Nettelmann 2010; Kang et al. 2011, 2013; French et al. 2012; Militzer 2013), which 
includes transport properties of the characteristic materials. In fact, atomic diffusion, 
as a basic physical process, plays a signihcant role in modeling the interior evolution of 
astrophysical objects (Michaud et al. 2007). Therefore, the physical properties of helium 
under extreme conditions is of continuous and strong interest owing to the aforementioned 
wealth applications to astrophysics (Fortov 2011). 

Despite the simplicity of atomic structure, bulk helium exhibits complex behavior, 
especially under astrophysical conditions because the states cover from atomic, partially 
ionized to fully ionized helium (Kietzmann et al. 2007; Khairallah & Militzer 2008; Stixrude 
& Jeanloz 2008; McMahon et al. 2012). In past few decades, great progresses have been 
achieved to explore the physics of helium at extreme states experimentally, such as diamond 
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anvil cell and shock wave (Celliers et al. 2010; Soubiran et al. 2012). However, only narrow 
bands of temperature and pressure parameters attainable by laboratory measurements 
today. Thus the comprehensive understanding of helium at extreme states requires accurate 
theoretical calculations. Conventional approaches that are suitable for treating ideal plasma 
gas or condensed matter under ambient conditions usually fail to solve the issues in the 
warm dense regime (Graziani et al. 2014). Ab initio approaches based on density functional 
theory (DFT) provide a feasible and satisfactory tool to investigate the physical properties 
of matter under extreme conditions. They have been extensively used to obtain the 
thermodynamical, transport, and optical properties of matter in astrophysical area (French 
et al. 2012). Quantum Langevin molecular dynamics (QLMD), including non-adiabatic 
effects such as the electron-ion collision induced friction in molecular dynamics (Dai et al. 
2010a), extends the quantum molecular dynamics to the regime up to solar interior (Dai et 
al. 2010b, 2012, 2013). 

A common feature of these theoretical methods is that the ions are treated classically, 
that is, the quantum nature of ions motion is neglected. It should be noted that due to 
the low mass, helium motion is expected to be fairly sensitive to quantum effects, the 
so-called nuclear quantum effects (NQEs), as hydrogen exhibits (Kang et al. 2014). The 
well-known superfluidity of helium is dominated by quantum effects of helium atoms at 
extremely low temperatures (Ceperley 1995). At relatively high temperatures in planet’s 
core or outer layer of WDs, NQEs are usually considered to be negligible before. Note that 
with the increasing density of helium, the atomic distance decreases. When the nearest 
atomic distance is comparable to the thermal de Broglie wavelength, the nuclear quantum 
delocalization would play remarkable roles in structure and transport properties, e.g., NQEs 
induce complex transport behavior in warm dense hydrogen as reported in our previous 
work (Kang et al. 2014). Path-integral molecular dynamics is an effective tool to investigate 
the NQEs on static properties of extreme states (Marx & Parrinello 1996) and recently 
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improved centroid path-integral molecnlar dynamics (PIMD) can provide the well-dehned 
qnasiclassical approximation of real-time evolntion of qnantnm systems (Cao & Voth 1993; 
Kang et al. 2014). Therefore, centroid PIMD, combined with DPT, treats both electrons 
and ions qnantnm-mechanically. Using this method, more accnrate description of transport 
properties of matter at extreme states is expected to be obtained. 

In this stndy, we investigate NQEs on the strnctnral and transport properties of dense 
liqnid helinm at 5 g/cm^ and 10 g/cnU nsing ab initio centroid PIMD with an improved 
scheme, and compare the resnits with that from traditional hrst principles molecnlar 
dynamics (MD) calculations based on DFT with a classical treatment of nuclei to assess 
the importance of NQEs. 


2. Theoretical Methods 

Path-integral formulation of quantum statistical mechanics provides a theoretical and 
computational framework to study the quantum many-body system at hnite temperatures 
(Feynman 1972). In path-integral picture, the canonical partition function can be expressed 
as a conhgurational integral of Boltzmann weighted continuous paths. If exchange effects of 
particles are negligible, these closed paths can be discretized through Trotter approximation 
into P beads circularly connected via harmonic springs, so that a quantum system including 
N particles is isomorphic to a classical system consisting of N ring polymers, where each 
quantum particle is mapped onto a closed flexible polymer of P beads (Parrinello &: 
Rahman 1984; Marx &: Parrinello 1996). Static properties, e.g., radial distribution function 
(RDF), of the quantum system are routinely obtained through molecular dynamics sampling 
the conhguration space of its isomorphic classical system. Further, centroid molecular 
dynamics within the framework of path-integral improved recently is a good qnasiclassical 
approximation to quantum dynamics (Kang et al. 2014). In the adiabatic centroid PIMD 


scheme (Marx et al. 1999), the primitive path variables (s is the imaginary time.) 

is transformed to a set of normal mode variables which diagonalize the harmonic 

bead coupling. The equations of motion of normal mode variables without thermostats are 
given by 



( 1 ) 



^yi s'=i 



( 2 ) 


where up = y/PkpT, P is the number of beads in ring polymer, kp is Boltzmann constant. 


T is the target temperature of the simulated system, is the electronic energy functional 


in ab initio calculations. 

The centroid and non-centroid modes are mass-scale separated with adiabatic 


parameter 7 , i.e., the normal mode masses are 


= Mp , s = 2, • • • , P. 


( 3 ) 


where Mj is the physical nuclear mass. In this way, the nuclear quantum dynamics is 


obtained in the quasiclassical approximation through driving the centroid move in real-time 
in the centroid effective potential generated by all non-centroid modes. Here we propose 
a new scheme for sampling the force on the centroids during the molecular dynamics 
simulations to make the size of ring polymer is more close to the de Broglie wavelength 
as the instantaneous kinetic energy of ions is decreased. That is, for each ion ujp depend 
linearly on its instantaneous kinetic energy E{t) during simulations. 



( 4 ) 


9 


where t is the simulation time step. The details of determining the parameters a and b are 
given otherwhere (Kang et al. 2014). Here we show the results, i.e., a = 1 — (^^)^ and 
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b = where RgQ is the radius of gyration of ring polymer when E(t) = Eq, X is the de 

Broglie wavelength of the corresponding ion. As shown in Kang et ah (2014), this improved 
centroid PIMD approach can describe the quantum collision properly, and nuclear quantum 
dynamics in dense hydrogen has been investigated successfully. 


3. Results and Discussion 
3.1. Computational details 

For dense helium, we study the structural and transport properties at densities of 5 
g/cm^ and 10 g/cm^ above the melting point temperature. The knowledge of the structure 
and transport properties of dense liquid helium under these density and temperature 
conditions is essential to correctly describe the interior structure and evolution and cooling 
of astrophysical bodies, such as exoplanets and outer layer of white dwarfs. Exchange effects 
between ions are neglected in this study (Dyugaev 1990) because the atomic distance is still 
larger than the thermal de Broglie wavelength of ions under the conditions considered here. 

Our MD and PIMD simulations were performed using the modihed Quantum- 
ESPRESSO package (Giannozzi et al. 2009), where electrons are quantum-mechanically 
treated by hnite-temperature density functional theory (Mermin 1965). The interaction 
between valence electron and the ionic core is presented by norm-conserving pseudopotential. 
The generalized-gradient approximation for exchange-correlation potential (Perdew et 
al. 1996) is adopted. The plane-wave kinetic energy cutoff was set from 150 to 180 Ry 
according to different densities. The Brillouin zone is sampled using the Gamma point in 
dynamic simulations, while denser k-point grids of 4x4x4 Monkhorst-Pack scheme points 
are used to calculate electronic properties. The periodic supercell including 128 or 250 
atoms is adopted for 5 g/cm^ and 10 g/cm^, which can ensure convergence of both ionic 



and electronic properties with good accuracy. 


Within the framework of PIMD, the structural properties are calculated using the 
primitive scheme, while the real-time quantum dynamics of nuclei is obtained through 
the improved centroid PIMD. The self-diffusion coefficient is calculated from the slope 
of the centroid mean square displacement, while the shear viscosity is obtained from the 
autocorrelation function of the off-diagonal components of the stress tensor (Alfe & Gillan 
1998). Langevin thermostat is employed to overcome the nonergodic problem, which not 
only produces a canonical ensemble and compensates the calculated errors, but also gives an 
good method to include the non-adiabatic effects of electrons in warm and hot dense regime 
(Dai et ah 2010a,b). Note that Langevin thermostat is applied only to each noncentroid 
degree of freedom because large thermostats would disturb the dynamical properties of the 
centroids. The Trotter number was set to 16 after a convergence test. Time step of 7-9 a.u. 
is used in the MD calculations with 30000 steps after thermalization, while the smaller time 
steps of 0.3 a.u. with 5x10^ steps and a centroid adiabaticity parameter of 20 are adopted 
in the centroid PIMD simulations in order to decouple the centroid mode from other normal 
modes. 


3.2. Transport properties of dense liquid helium 

Firstly, we calculated the RDF of helium nuclei to reveal the structural difference 
introduced by NQEs. Within the framework of imaginary-time PIMD, the quantum 
expectation value of a static physical quantity is given in terms of averaging over the 
canonical ensemble of the isomorphic classical system. Therefore, the radial distribution 
functions of helium atoms are obtained by averaging over both the PIMD time steps and 
imaginary time slices. Figure 1 shows that NQEs play an important role in RDF at low 
temperatures with the density range of 5 g/cm^ and 10 g/cm^. In particular, the hrst peak 
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of RDF from PIMD simulations are largely reduced and broadened, which indicates the 
significant nuclear quantum delocalization. When the temperature is increased, the nuclear 
quantum delocalization becomes weaker and cannot be observed ultimately at 10000 K for 5 
g/cm^ helium and 12000 K for 10 g/cm^ helium. Meanwhile, RDF shows more pronounced 
nuclear quantum effect with increasing density, even though the temperature (7000 K) at 
10 g/cm^ is higher than that (4000 K) at 5 g/cm^. 

Of our particular interest is NQEs on transport properties of dense liquid helium 
because atomic diffusion is the basic physical process that determine the interior structure 
and evolution of astrophysical objects (Michaud et ah 2007). The results of the self-diffusion 
coefficient and shear viscosity at 5 g/cm^ with the temperature range of 4000-10000 K 
are shown in Figure 2. When taking into account the quantum features of the collision 
processes between ions properly, the self-diffusion coefficients obtained by centroid PIMD 
simulation are substantially larger than the classical simulation value over the whole 
temperature range of 4000-10000 K. The difference is 37% at 4000 K, 22% at 6000 K, 
17% at 8000 K and 10% at 10000 K. It is interesting that even though the RDF from MD 
and PIMD calculations are almost identical when the temperature is increased to 10000 
K, the ionic diffusion still exhibits distinct difference with the inclusion of NQEs. In fact, 
the essence of atomic diffusion is atomic scattering from the viewpoint of collision physics. 
The large-angle scattering between ions is dominant under the high density conditions 
considered here. The pronounced quantum nuclear scattering features introduced by 
centroid PIMD simulations leads to a smaller large-angle scattering cross section between 
ions than the classical treatment, thereby increasing the mean free path of ions. Thus 
the classical-particle treatment of protons substantially underestimates the ionic diffusion. 
Similarly, the self-diffusion coefficients at the temperature range of 7000-12000 K from 
quantum simulations for helium at 10 g/cm^, as shown in Figure 3, are 23 %, 21 %, 12 % 
and 6 % higher than that from classical simulations. 
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The shear viscosities calculated from the Green-Kubo relation (Alfe & Gillan 1998) 
at the density of 5 g/cm^ and 10 g/cm^ are also shown in Figure 2 and 3, respectively. 
We can see that the viscosities from centroid PIMD are smaller than those from MD and 
the difference between them becomes smaller with increasing temperatures. In addition 
to obtaining the viscosity from molecular dynamics simulations, it can be alternatively 
deduced from the Stokes-Einstein (SE) relation with slip or stick boundary conditions 
(Hansen & McDonald 2011). The SE relation, which is exact for the Brownian particles, is 
not valid in strongly coupled regime. In order to examine the SE relation, we calculated 
the parameter Dr^a/ksT (a is the effective atomic diameter and dehned by the position 
of the hrst peak of RDF) and the results are shown in Figure 2 and 3 for the density of 
5 g/cm^ and 10 g/cm^, respectively. The parameter Dria/ksT both from MD and PIMD 
calculations deviates from the SE relation value l/ 27 r at low temperatures; furthermore, 
NQEs introduce pronounced difference between the classical and quantum simulations. In 
fact, the state of helium in planet’s core or cool atmosphere of WDs is in moderately or 
strongly coupled regime and the motion of ions is quite different from Brownian motion. 

In order to understand the atomic diffusion behavior with NQEs microscopically, we 
explore the potential surface of a randomly selected atom along its simulation trajectory. 
Given the force in simulations is calculated by F' = —VR, the potential energy along the 
trajectory of a atom can be obtained hj V = Vq + J F ■ dr, where F is the force acting 
on atoms, V is the potential energy, Vq is the reference point of potential energy. The 
potential energy of a randomly selected helium atom along its simulation trajectory from 
MD and PIMD simulations as well as the distribution of the potential energy is presented in 
Figure 4. We can see that in PIMD simulations the potential energy of the imaginary-time 
beads fluctuates steeply, which presents the effects of quantum fluctuations on dynamics 
of the centroid. The potential surface of the centroid exhibits more smoothly than that 
of classical particle, which is evidenced by the fact that the distribution of the potential 
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energy of centroid covers smaller energy range than MD simulation and the high frequency 
component of the potential surface of the centroid is higher than MD results by about 100 
THz for both 5 g/cm^ and 10 g/cm^, as shown in Figure 4(c) and (f). It indicates that when 
considering NQEs the atoms are more prone to turn over the energy barrier than classical 
treatment, thereby enhancing the atomic diffusion. In addition, it is shown in Figure 4 that 
when the density is increased from 5 g/cm^ to 10 g/cm^, the potential surface exhibits more 
barriers and wells in the same range of time (100 fs) both in MD and PIMD simulations. It 
unambiguously indicates that helium atoms collide more frequently at higher density at the 
same temperature so that the atomic diffusion is lowered at high densities. 

Nuclear motion quantum nature not only affect the atomic diffusion behavior, but also 
induce the electronic redistribution (Cannuccia & Marini 2011), thus introducing more 
localized electrons (Kang et ah 2013, 2014). It can be evidenced by the distribution of 
charge density shown in Figure 5. Note that with the inclusion of NQEs, the charge density 
has the larger maxima and smaller minima compared with the classical simulations, which 
indicates that electrons surrounding ions distribute more localized with quantum nuclei 
than classical particle treatment. For metallic state of matter, e.g., dense hydrogen, this 
electronic localization can signihcantly lower the dc electrical and thermal conductivities 
(Kang et ah 2014). For dense liquid helium, NQEs induce a strong impact on the electronic 
density of states (DOS). We can see from Figure 6 that the energy bands, especially 
conduction bands are obviously broadened when NQEs are considered via PIMD. Since 
helium exhibits a large electronic energy gap under the conditions considered here, it is still 
in insulator state. Therefore, NQEs on electronic transport properties are not discussed 
here. Except for NQEs on DOS, density effect has also influenced the prohle of DOS 
remarkably. As shown in Figure 6, when the density of helium is increased from 5 g/cm^ to 
10 g/cm^, the energy gap is decreased, thus the K-edge photo-absorption spectra will have 
a red shift. More importantly, we note that the DOS of the Is state is extended by about 
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30 eV when helium is compressed from 5 g/cm^ to 10 g/cm^, which signihcantly affects the 
interaction potential between helium nuclei and make the potential surface fluctuates more 
frequently at high density (shown in Figure 4). Therefore, both NQEs and density effect 
contributes to the complex transport properties of atoms. 


4. Conclusion 

In conclusion, NQEs on transport properties of dense liquid helium under the conditions 
of the planet’s core and cool atmosphere of WDs have been investigated using the improved 
centroid PIMD simulations combining with DFT. The results show that with the inclusion 
of NQEs, the self-diffusion is largely higher while the shear viscosity is notably lower than 
the results of without the inclusion of NQEs due to the lower collision cross sections even 
when the NQEs have little effects on the static structures. Meanwhile, the relation between 
diffusion and viscosity is deviate from the SE relation valid for Brownian particles, thus 
the SE relation is not valid in strongly coupled regime. The potential surface of helium 
atom along the simulation trajectory is quite different between MD and PIMD simulations. 
We have shown that the quantum nuclear character induces complex behaviors for ionic 
transport properties of dense liquid helium. Therefore, in order to construct more reasonable 
structure and evolution model for the planets and WDs, NQEs must be reconsidered when 
calculating the transport properties at certain temperature and density conditions. 

This work is supported by the National Basic Research Program of China (973 
Program) under Grant No. 2013CB922203, the National Natural Science Foundation of 
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Fig. 1.— Radial distribution functions of helium nuclei from PIMD (solid lines) and MD 
(dashed lines) simulations with different temperatures at 5 g/cm^ (a) and 10 g/cro? (b). 
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T(K) 

Fig. 2.— Temperature dependence of the self-diffusion coefficient and shear viscosity at 5 
g/cm^ with temperatures of 4000 K, 6000K, 8000 K and 10000 K. obtained from centroid 
PIMD and MD simulations. The shear viscosities obtained from the Stokes-Einstein (SE) 
relation are also presented. The effective atomic diameter in SE relation is dehned by the 
position of the hrst peak of radial distribution function. 
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T(K) 


Fig. 3.— Temperature dependence of the self-diffusion coefficient and shear viscosity at 10 
g/cm^ with temperatures of 7000 K, 8000K, 10000 K and 12000 K. obtained from centroid 
PIMD and MD simulations. The shear viscosities obtained from the Stokes-Einstein (SE) 
relation are also presented. The effective atomic diameter in SE relation is dehned by the 
position of the hrst peak of radial distribution function. 
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Fig. 4.— The potential energy of a randomly selected helium atom along its simulation 
trajectory. The results of classical simulation (dashed lines), the centroid of PIMD (solid 
lines) and the randomly-selected bead in ring polymer (dotted lines) are displayed in two 
conditions of (5 g/cm^, 8000 K) and (10 g/cm^, 8000 K). Middle panels (b) and (e) are 
the distribution of the potential energy. Note that the peaks of distributions are shifted to 
the same value for comparisons. Right panels (c) and (f) are the fourier transform of the 
potential surface. 
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Fig. 5.— The cut-plane of 3D distribution of charge density of randomly selected conhg- 
urations from MD (a, c) and PIMD (b, d) simulations. For upper panel, the density is 5 
g/cm^ and the temperature is 4000 K. For lower panel, the density is 10 g/crv? and the 
temperature is 7000 K. In PIMD calculations, the conhguration is randomly selected from 
imaginary-time slices in simulation time steps. 




Fig. 6.— The electronic density of states (DOS) calculated from quantum (solid lines) and 
classical simulations (dashed lines) at 5 g/cm^(a) and 10 g/cm^(b). The results from PIMD 
are averaged over not only the simulation time step but also the imaginary-time slices in 
each simulation time step. 








